Observation of String Breaking in QCD 
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We numerically investigate the transition of the static quark-antiquark string into a static-light 
meson- ant imeson system. Improving noise reduction techniques, we are able to resolve the signature 
of string breaking dynamics for n/ = 2 lattice QCD at zero temperature. This result can be related 
to properties of quarkonium systems. We also study short-distance interactions between two static- 
light mesons. 
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I. INTRODUCTION 

Sea quarks are an important ingredient of strong in- 
teraction dynamics. In the framework of quantum chro- 
modynamics, however, quantitative calculations of their 
effects on hadron phenomenology have proven to be no- 
toriously difficult, unless one resorts to approximations 
based on additional model assumptions. Nevertheless, 
the ab initio approach of lattice gauge theory towards 
the sea quark problem has shown steady progress over 
the past decade: recently the ?/-problem has been tackled 
successfully on the lattice [1-3] where sea quarks induce 
the axial anomaly in the sense of the Witten-Veneziano 
mechanism [4, 5]. 

Another example is the strong decay of hadrons 
through light quark-antiquark pair creation, for instance 
the transition from a colour string configuration between 
two static colour sources, QQ, into a pair of static-light 
mesons, BB. This colour string breaking, which we ad- 
dress in this paper, is expected to occur as soon as the 
colour source-sink separation, r, exceeds a certain thresh- 
old value, r c > 1 frn. 

In lattice simulations this behaviour has been investi- 
gated in four dimensional QCD at zero temperature T 
with sea quarks [6-12] as well as in QCD3 [14]. However, 
these studies lacked compelling evidence of string break- 
ing 1 . This failure is due to problems like: (i) String 
breaking investigations only make sense in a full QCD 
setting with large ensemble sizes, (ii) String breaking 
occurs at distances beyond 1 fm, a regime with a poor 
signal-to-noise ratio. (Hi) The poor overlap of the QQ 
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1 The T > situation appears to be more favourable [13]. 



creation operator with the large-distance BB ground 
state. 

This last problem necessitates to resolve the signal at 
huge Euclidean times t, unless one bases the investigation 
on a 2 x 2 correlation matrix, whose additional elements 
include the insertion of light quark propagators into the 
standard Wilson loop [7, 10-12]. Such quark insertions 
require propagators from any source to any sink position 
("all-to-all propagators"), in order to enable the exploita- 
tion of translational invariance for error reduction (self 
averaging). 

For QCD with n f mass-degenerate sea quark flavours 
this correlation matrix takes the form, 



C(t) = 



( C QQ (t) Cqb (t) 
\C BQ (t) Cbb (t) 

( 



\ 



(1) 



where the straight lines denote gauge transporters and 
the wiggly lines represent light quark propagators 2 . We 
refer to the difference between the physical eigenstates 
and the QQ and BB basis as "mixing". Such mixing 
should manifest itself "explicitly", by non- vanishing off- 
diagonal matrix elements, relative to the diagonal matrix 
elements, and "implicitly" . The latter refers either to 
the Wilson loop Cqq{€) decaying into the mass of the 
(dominantly) BB state for r > r c or to a decay of Cb b if) 
towards the QQ mass for r < r c , as t — > 00. Implicit 
mixing is much harder to detect than explicit mixing. 

In the quenched approximation baryon and anti- 
baryon numbers are separately conserved and the QQ 
and BB sectors are mutually orthogonal. By definition, 



Details of this expression will be discussed in Sec. II below. 
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rif vanishes and there will be no mixing. This does not 
mean however, that the matrix elements accompanying 
the explicit nj- and ^/nj-factors are necessarily zero. 

So far string breaking has been verified in the following 
cases: in SU(2) gauge theory with a fundamental scalar 
field in three dimensions [15] and in four dimensions [16- 
18] as well as for the SU(2) potential between adjoint 
sources (screened by the gluons) in three dimensions [19- 
22] and in four dimensions [23, 24]. However, in only one 
of these studies [22] , and in a recent simulation of the 3d 
Z 2 -Higgs model [25], implicit string breaking has been 
convincingly demonstrated. 

Let us recall the signature of string breaking: without 
mixing, the QQ and the BB are QCD eigenstates and 
will undergo a plain level crossing (with minimal energy 
gap AE C = 0), at a certain critical distance r c . Sim- 
ulations with nf = show exactly this behaviour. In 
contrast, with qq creation/annihilation switched on, the 
Fock states will undergo sizeable mixing in the neigh- 
bourhood of r c . The minimal energy gap AE C between 
the two eigenstates will grow with the spatial width of 
the mixing region. 

The transition rate between QQ and BB 
states is given by the (normalized) time deriva- 
tive of the off-diagonal matrix element, g — 
[dC Q B(t)/dt} t=0 [CBB(0)C QQ (0)}-^ 2 . From a 'string 
picture as well as from strong coupling arguments one 
would expect a more pronounced mixing in larger space- 
time dimensions d. Therefore, the size of the energy gap 
within the string breaking region should increase as one 
goes fro m d = 3 to d = 4. In the large -ZV C limit we find, 
g oc y/nf/N c , for the potential between fundamental 
sources, screened by n/ flavours of fundamental scalar 
or quark fields. For the breaking of the adjoint potential 
this translates into, g oc 1/N C . 

The figures of Refs. [20, 21] for the adjoint string in 3d 
SU(2) gauge theory suggest the following upper limits 
for the size of the energy gap, expressed in units of the 
string breaking distance: AE c r c < 0.45 and AE c r c < 
0.75, respectively. In 4d SU(2) gauge theory one finds, 
AE c r c < 0.6 [23], whereas Refs. [15, 18] show that the 
fundamental 3d SU(2) string, screened by a scalar field, 
satisfies AE c r c < 0.35 and AE c r c < 0.65, respectively. 
In all these cases either the spatial resolution of the string 
breaking region was too coarse or the statistical errors 
were too large to allow for the determination of a lower 
bound. 

Based on the qualitative n/, N c and ^-dependencies 
discussed above, we expect the n/ = 2 QCD energy gap 
to be somewhat bigger than the gaps quoted for the toy 
model studies. From this reasoning we would aim at an 
error of AE c r Cl smaller than 0.1. To meet this constraint, 
we require a distance resolution in the string breaking 
region of Ar < 0.1/(or c ) w 0.02r c w 0.025 fm, where 
a denotes the string tension. We will find, AE c r c — 
0.33(5). 

In order to achieve the required precision, we apply a 
fourfold arsenal of critical improvements, within the 2x2 



correlation matrix setting: 

Ground state overlaps. It is essential to achieve a 
large overlap between the trial wavefunctions and 
the respective physical ground states. This en- 
hances the signal since it will decay less rapidly 
with Euclidean time. Moreover, the large t asymp- 
totics will be reached at smaller temporal distances, 
further reducing the noise/signal ratio. To this end 
we employ combinations of APE and Wuppertal 
smearing techniques (see Sec. IIIB). 

Wilson loops. The Wilson loop signal Cqq^) can be 
further enhanced by using an improved fat link 
static action. In this way, the relative errors of the 
Wilson loop data are reduced by factors of about 
five (see Sec. IIIB). 

Quark propagators. The generalized Wilson loops of 
Eq. (1) require the computation of all-to-all light 
quark propagators if we wish to fully exploit self 
averaging. Direct inversion of the Wilson Dirac 
matrix M would be computationally prohibitive. 
Therefore, we approximate M~ Y by the lowest lying 
eigenvectors of 75M using the truncated eigenmodc 
approach (TEA) [26] , together with a stochastic es- 
timation (SET, see e.g. the review [27]) in the or- 
thogonal subspace. Moreover, we apply a "hopping 
parameter acceleration" (HPA) for variance reduc- 
tion. This further reduces the errors of the dis- 
connected contribution to Cb b by factors of about 
three (see Sec. IIIC). 

Distance resolution. In order to avoid finite size ef- 
fects and to achieve a fine distance resolution of 
the string breaking region, we employ a large set of 
off-axis distances (see Sec. Ill A). 

With these methods we are able to demonstrate com- 
pelling evidence, both for explicit mixing and for string 
breaking in full QCD, as well as for implicit mixing within 
CsB{t) for r < r c . We find the breaking of the quark an- 
tiquark string to occur at a distance r c w 15 a « 2.5 r w 
1.25 fm, in units of r « 0.5 fm [28, 29]. 

Note that our study should be viewed as exploratory 
since we restrict ourselves to one value of the sea quark 
mass (slightly below the physical strange quark mass), 
at one lattice spacing. For details on our simulation pa- 
rameters, see Sec. Ill A. 

The paper is organized as follows: in Sec. II we dis- 
cuss the mixing problem in detail. The notation used 
within Eq. (1) will be defined. We describe the com- 
bined application of TEA, SET and HPA for the calcu- 
lation of all-to-all propagators in Sec. III. In Sec. IV we 
discuss theoretical expectations for the individual matrix 
elements and check these against our data. In Sec. V, we 
present and interpret our main result, string breaking in 
QCD, as well as the short-distance interactions between 
two static-light mesons with isospin I — and 1=1. 
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We comment on the phenomcnological implications in 
Sec. VI. 

In view of the length of this paper, we kept the sections 
as self-contained as possible. For instance, the reader 
who is less interested in the technical aspects of the study 
can safely skip Sec. Ill altogether and concentrate on 
Sees. V and VI. 



II. THE MIXING PROBLEM 

Let us consider a system with a heavy quark Q and 
a heavy antiquark Q in the static approximation, with 
separation r = |R|/a, where a denotes the lattice spacing 
and R is an integer valued three-vector. We restrict our 
discussion to the E+ ground state of the static system, 
with cylindrical symmetry. 

Without sea quarks, the energy of this static-static sys- 
tem will linearly diverge with r [30-32] as r — > oo. In the 
presence of sea quarks, however, there will be some crit- 
ical "string breaking" distance r c : when r exceeds r c , 
the mass of a system containing two static-light mesons, 
which we shall call B and B, separated by r, will be- 
come energetically favoured. The static QQ potential 
will exhibit screening and saturate towards about twice 
the mass of the B meson. 

A full investigation of this phenomenon requires the 
study of the Green functions that correspond to the prop- 
agation of the QQ and BB systems as well as of the 
transition element between these two states. We start by 
defining our notations and discussing the symmetries of 
the problem, before we display the relevant Green func- 
tions. 



A. Definitions and representations 

The Euclidean Dirac equation in the static limit, 

(D AlA + m Q )\Q) = 0, (2) 
yields the propagators of static quark and antiquark: 

Q y Q x = (y\Q)(Qhi\x) (3) 

Q f y Ql = (y\Q)(Q\i4\x) (4) 

where y 4 > x 4 and, 

1±74 



(5) 



are projectors onto the upper and lower two Dirac compo- 
nents. t/x(i/4, £4) e SU(3) denotes a lattice discretization 
of the Schwinger line, connecting (x, X4) with (x, y 4 ): 



T denotes the time ordering operator. We use the conven- 
tion (x\Q) = Q x and (Q\x) = Ql, i.e. Q x = (Q|7 4 |a;). 

ttiq (a) in Eqs. (3) and (5) above is the heavy quark 
mass in a lattice scheme 3 . We define the light quark 
Dirac operator, 

M = (\q)(q\ l4 y 1 = D^+m. (7) 

We use Wilson fermions throughout the paper: 

4 

M X y — S X y — K ^ ' [(1 — "f^Ux^Sx + aft^y 

x — a£i,y ]• (8) 

As usual, the quark fields have been rescaled by factors 
\^2k where n = (8 + 2ma)~ 1 in the free field case and in 
general, n(m = 0) = k c > 1/8. U x ^ denotes an SU(3) 
gauge field and U x _ p = £^_ aAijU , Aft = 75M75. 

We define a gauge transporter Vt(y,x) <G SU(3), con- 
necting the point (x, t) with (y,i). This is to be taken 
local in time and rotationally symmetric about the short- 
est connection. The properties under local gauge trans- 
formations tt x E SU(3) are: 



Qx 1 > QxQxj Qx 1 > Qx^xi 
V t (y,x) 1 ^ fi(y,t)T4(y,x)fi| Xit)) 

which implies that the combination, 



Qt(y,x) = Q 



o,*)^r^(y> x )Q(x,t), 



(9) 
(10) 



(ii) 



is a colour singlet. 

The spins of Q and Q can either couple symmetrically 
or anti-symmetrically. The first situation is represented 
by 7 • r/r in Eq. (11) (total spin 5=1), the latter choice 
corresponds to the replacement, 7 • r/r 75, within 
Eq. (11) (5 = 0), see also Ref. [35]. 

The relevant symmetry group is not 0(3) ® C but its 
cylindrical subgroup Dooh. On the lattice this reduces 
to D^h- Nevertheless, we will use the continuum expres- 
sions, as the "latticization" is straight forward in this 
case [36, 37]. The irreducible representations of D^h 
are conventionally labelled by the spin along the axis, A, 
where E, II, A refer to A = 0,1,2, respectively, with a 
subscript 77 = g for CP = + (gerade, even) or 77 = u for 
CV = — (ungerade, odd) transformation properties. Par- 
ity V or charge C alone are not "good" quantum numbers. 



jig J dt A 4 (x,i)j . 



U x {y 4 ,x 4 ) ~ T cxp \ig I dt Ai(x,t) \ . (6) 

/ x& 



3 rriQ contains a power term in the inverse lattice spacing o _1 , 

<5m k ai/o H , where «i = g /{Air) is the strong coupling 

parameter in the lattice scheme. This term (that diverges in the 
continuum limit a — > 0) cancels against a similar contribution 
from t/ x (t2,ti). We shall also refer to this contribution as the 
"self-energy" associated with the static propagator. Note that 
factorising rug (a) into pole mass and self-energy introduces a 
renormalon ambiguity, see e.g. [33, 34]. 
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The £ representations carry, in addition to the r\ quan- 
tum number, an 1Z parity with respect to reflections on a 
plane that includes the two endpoints. This results in an 
additional ± superscript. The symmetric spin combina- 
tion Eq. (11), when combined with a symmetric gluonic 
string Vt, lies within the £+ ground state representation 
while the antisymmetric spin-combination, 7 • r/r 1— ► 75 
corresponds to S~. These two representations yield de- 
generate energy levels, since both are calculated from one 
and the same Wilson loop 4 . 

Once mass corrections are added to the static limit, 
the full 0(3) £g> C symmetry becomes restored. The 
representation is contained within the J PC — 
++ , 1 , 2 ++ , • • • sectors of this bigger symmetry group 
while T,~ corresponds to J PC = h , 1 H Within 
the two-quark sector, the 1 and h states form the re- 
spective (mass-degenerate) ground states since the other 
quantum numbers require angular momentum L > or 
non-trivial gluonic excitations. 

B. The elements of the correlation matrix 

We consider nj mass-degenerate flavours of light 
quarks q l , i = 1,. ..,n/. Let Bi be the S-meson with 

4 The labelling is somewhat different if we start from scalar rather 
than from fcrmionic static sources. In this case the ground state 
E^" potential is not accompanied by any other mass-degenerate 
states and T,^ would label a non-trivial gluonic excitation. 

5 One can also work out the correlation matrix elements, starting 
from two vector states. Another possibility would be to probe 



light quark flavour i. For simplicity we label the Q r Qo 
string creation operator as Q and the Bi r Bi_o operators 
as Bi — Bu. We suppress the distance r for ease of nota- 
tion. We define the (unnormalized) states, 



|Q) = Q|0), |B j )=iB i |0). (12) 



In what follows, \Q) will always denote this QQ state 
and should not be confused with the static quark spinor 
of the same name within Eqs. (2) - (5). The lightest 

static-light meson has light quark J p = i , see e.g. [38] 
and references therein. Combining this with the heavy 
quark spin leads to mass-degenerate pseudoscalar and 
vector states. Two of these pseudoscalars/vectors com- 
bined have CP = + and fall into the £+ representation 
(in the vector case the spins have to be anti-aligned ac- 
cordingly, to yield J z = 0). For concreteness we shall 
choose 5 , Bi — Q"f§q l . We are now in the position to dis- 
play the three Green functions that are relevant to our 
problem. For the time evolution of the QQ state this 
reads: 

the sector with an antisymmetric combination of vector and 
pseudoscalar B states. As it should be, all these starting points 
yield identical Green functions (with the exception of r = 0), in 
the infinite quark mass limit. 



{q\t\q) v = (0 



T 



= 2e~ 2m «*(tr 



Q{T,try-V t (r, 0)Q (0it) ] f Q {rfi) ^V Q (r, 0)Q (0 , 0) 



_ e —2m,Qt 



{v*(T,0)U T (t,0)V (T,0)U*(t,0)}) u 



(13) 



where (0)u denotes the expectation value of O over 
gauge configurations. 

The trace above is over colour only. Hence the normal- 
ization of the Wilson loop is, (W(r,0)) u = 3. T = e~ aU 
denotes the transfer operator and r = t/a the Euclidean 
time difference in lattice units. The factor two originates 
from tr{P + 7jP_7j} = S^tr Pf = 2%. 

Note that in the first step of the derivation, a mi- 
nus sign from the commutator [Q^iC^]^ = —Q2liQi 
is cancelled when permuting one quark held through the 



remaining three others, prior to the Wick contraction. 

Expectation values and correct pre-factors are under- 
stood to be implicit in the pictorial representation of the 
correlators. In this case the time direction is assigned 
to be vertical, the spatial separation to be horizontal. 
Light quark propagators will be represented as wiggly 
lines, static quark propagators and gauge transporters 
are shown as straight lines. 

Next we consider the transition clement between QQ 
and BiBi states: 



■5 



(B i \T T \Q) u = (o|r{Q (0>t) 75^ 0>t) 4 it) 75Q(r, t )Q M )^^o(r,0)Q (0 ,o)}|o) 
= e - 2 "^(Tr {p_^M ( -^. (rit) £/ P (t,0)Vo(r,0)i;t( t ,o)})^ 

= (Q\T T \B l ) u , 



e -2m Q t 














WW 



(14) 



where the trace above is over colour and Dirac indices 
and we have made use of the relations, 75P+ = P-75, 
7i75 = ~757i, H p + = P-li and P - = p -- 



Finally, the BiBi sector reads: 



{B i \T T \B j ) u = (0 |7 1 {g (0it )75g(o,t)9(r,t)75Q(r,t)'9( r ,o)759( r!0) g( , 0) 75Q(o,o)}| 

= e - 2mQt % (Tr {P-M^ ty{rO) U T (t,0)}Tr {P + M^ oy{ot) Ul(t,0)} 

-e~ 2mQt {^-M ( - o ^ );(ri ^ r (i,0)P + M ( - 1 0);(Oi0) f/t (t)0 )}) [/ 
= e - 2m «* I 6, 



(15) 



Again, the traces are over colour and Dirac indices. In 
what follows, we will also refer to the flavour singlet sec- 
tor (which is the one relevant for the string breaking 
problem) as the 1 = sector while we label flavour non- 
singlet states as / = 1. In the isospin 1 = sector the 
connected diagram always contributes while the discon- 
nected diagram only contributes for i = j. In contrast, 
within the 1=1 sector there is no connected diagram 
but only the disconnected contribution 6 . 



C. Reduction to a 2 x 2 matrix 

We consider the scenario with rif > 1 mass-degenerate 
quark flavours. By summing over the flavour indices i 
and j in the above equations, the correlation matrix can 



If we allow the heavy quarks to move, by adding a kinetic term 
within a Born Oppcnheimcr approximation, then the 7 = sector 
can be related to transitions between vector bottomonium (or 
0++,2++, . . . bottomonia) into a pair of B and B pscudoscalar 
mesons or into B* and B vector mesons (which are mutually 
degenerate in mass in the static limit), with relative angular 
momentum L chosen appropriately. The simplest example for 
this sort of process is T(45) — ► BB with final state L = 1. 
Within the 7 = 1 sector one can write down a similar mixing 
problem. In this case, the other state would be a QQ plus an 
7 = 1 meson, like the n. The simplest such transition is BB —* 
Vb + 7T- 



effectively be reduced to a 2 x 2 problem: one can easily 
orthogonalize the BiBi meson-meson states; for given r 
and t all correlators, (Bi\T T \Bj)u, only involve one or 
two (for i = j) generalized Wilson loops [which are dis- 
played on the right hand side of Eq. (15)]. For rif = 2 
one can define for instance, 



\B) = ±={\B{) + \B 2 )), 
\B a ) = -J=(|Bi)-|B 2 ». 
Obviously \B a ) decouples from the other states: 



(Q\T T \B a )u = (B\T T \B a )u 
(B a \T T \B a )u = e- 2m « l 



0, 



(16) 
(17) 

(18) 
(19) 



This pattern easily generalizes to «/ > 2: as soon 
as two or more indices are antisymmetrized, the over- 
lap with the QQ state vanishes. Only the completely 
symmetric state has a non-trivial mixing (we write the 
formulas for general nj and \B) = -j== Ym=i 



(Q\r\B)u = ^WJe- 
(B\T T \B)u = e - 2m «* 




n f 



(20) 
(21) 



Hence, we have reduced the mixing problem for a gen- 
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eral nf to a 2 x 2 correlation matrix with elements, 

C af ,(t) = (a\T t / a \0) u , a,/3e{Q,B}. (22) 
This leads us to the form already anticipated in Eq. (1), 



C{t) = e 



_ —ImQt 



( 



V 



\ 



-n f 



(23) 

with the pictorial representations as defined in Eqs. (13) 
" (15). 

Note that we have some freedom to change the nor- 
malizations, without affecting the mass spectrum: 



Cqb >— > ab*CQB, Cbq i-> a*bCsQ, 



(24) 



C QQ ^ \a\ 2 C QQ , C BB >-> \b\ 2 C BB . (25) 

The phase of the Cqb = C* B q element is irrelevant and 
we have employed one of the two possible real choices. 
Cb b = Cg B — C B ° B consists of the following disconnected 
and connected contributions, 



Cdis 
BB 



-2'mqt 



rcon 
BB 



n f e 



-2mot 



• (26) 



In what follows we will set toq = 0, corresponding 
to a shift in all energy eigenvalues. Differences between 
two energy levels, such as between the mass of the QQ 
system and twice the static-light mass, do not depend 
on toq and have a well defined continuum limit. Wc 
remark, however, that the levels themselves become cut- 
off independent only in the framework of effective field 
theories. In this case, toq (a) is required to cancel the 
static self energy divergence and only the sum of quark 
masses and the potential is a "physical" quantity. 



III. MEASUREMENT TECHNIQUES 

We discuss the run parameters and the geometrical 
setup of our simulation, before wc elaborate on the noise 
reduction and all-to-all propagator techniques that we 
apply. We conclude by introducing the notations that 
we will use in the interpretation of our numerical data. 



A. Simulation set-up 

We base our simulations on the 24 3 x 40 T%L config- 
urations [39] with Wilson fermions at k = 0.1575 and 
[3 = 5.6. This translates into r = 6.009(53)a, cor- 
responding to a -1 w 2.37 GeV or a w 0.083 fm from 
ro = 0.5 fm. The value of tq differs somewhat from the 
earlier result, r = 5.892(27)a [6], that we obtained with- 
out accounting for mixing effects. While our results on 
Cqq(^) are much more precise than in this earlier study, 



after diagonalizing the mixing matrix, the final errors of 
the ground state energy level at r < tq increase. At larger 
r, however, we achieve unprecedented precision. One ob- 
tains [39] m^a = 0.276(5) and m^/my = 0.704(5), which 
means that the sea quark mass is slightly smaller than 
that of the physical strange quark. 

In order to stay clear of finite size effects, in partic- 
ular within the BB sector, it is advisable to place the 
colour sources off- axis. An on- axis string breaking study, 
in which r = n(l, 0, 0)a, n integer, would require a spatial 
lattice extent La > 2r c , r c being the string breaking dis- 
tance. Off-axis separations allow for a relaxation of the 
above condition to La > 2/^/3r c for the spatial diago- 
nal, r = n(l, 1, l)a, and La > 2j\f2r c along the planar 
diagonal, r = n(l, 1, 0). 

We have performed measurements on the following set 
of geometries: 



n(l,0,0)a,n < 11, 

n(l, l,0)a,n < 11, 

n(l, 1, l)a, n < 11, 

(10, 10, n)a,ne {2,3,4,5,6} 

(10,8,7)a, 

(10,9,7)a, 



as well as for r = 0. The distance resolution is en- 
hanced around r c w 15 a, to 10 points inside the range, 
14.14 w 10V2 < r/a < 15.59 w 9V3. In order to in- 
crease statistics we average over equivalent permutations 
and reflections of the axes. We do not find any directional 
dependence, even for the (11,0,0) and (11,11,0) points, 
and hence there is no sign of finite size problems close to 
r c . In the neighbourhood of r c , the largest component 
that we employ is r, = 10a. 

With r c w 15 a and L = 24 we obtain, V2i/(2r c ) w 
1.13 > 1 and V3i/(2r c ) w 1.39 > 1. Therefore, the T X L 
physical lattice extent La s=s 2.0 fm is sufficiently large for 
our purpose. We remark that the pion correlation length 
also fits well into the spatial lattice extent, La > 6m" 1 . 

We extract the elements of C(t) [Eq. (23)] that in- 
volve light quark propagators from a set of 20 thermal- 
ized gauge configurations {Ui}, i — 1, . . . , 20, separated 
by 125 Hybrid Monte Carlo trajectories. Earlier stud- 
ies [6, 40, 41] have established that these configurations 
are effectively independent of each other. 

The standard Wilson loop Cqq{€) is determined on a 
larger ensemble of 184 configurations, separated by 25 
trajectories. We also wish to eliminate possible autocor- 
relations in this case. Moreover, we attempt to consis- 
tently take account of correlations between different ma- 
trix elements (that have been determined on one and the 
same set of configurations). To this end, the 184 config- 
urations are averaged into 20 bins that are mapped onto 
the ensemble {Ui}. Each bin i contains the five configu- 
rations that are closest in Monte Carlo time to the above 
mentioned 20 configurations as well as an additional four 
to five configurations from within another region of the 
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time series. As it turned out, the limiting factor of our 
statistical resolution is the accuracy of the Wilson loop 
data and hence little can be gained from increasing our 
sample size for Cq_b(£) and CBs(i) beyond 20 configura- 
tions. 



B. Signal enhancement techniques 

We are interested in the exponential decay of the ele- 
ments of the correlation matrix C(t) at large Euclidean 
times. Statistically significant results cannot be achieved 
unless the asymptotic behaviour can already be extracted 
at moderate time separations. To this end, we employ 
smearing techniques that enhance the overlap of the op- 
erators used in the creation of particular states with the 
corresponding physical ground states, without affecting 
the eigenvalues. Furthermore, the noise/signal ratio has 
to be controlled. In pure gauge theories extended oper- 
ators can be constructed, with reduced variance, retain- 
ing identical expectation values [42, 43]. Unfortunately, 
these techniques, which exploit the locality of the gauge 
action in space-time, are not applicable when including 
sea quarks that (after integration) induce non-localities. 
Instead of reducing the variance, we enhance the signal 
by an appropriate choice of the lattice static quark ac- 
tion. 



1. Smearing 

We employ the iterative APE [44, 45] smearing pro- 
cedure for the spatial transporters V t that enter the cre- 
ation operator of the QQ states, Eq. (11): 



jAn+l) _ p 



5(7(3) 



(")t 
ai.j 



where i £ {1,2,3},.? e {±1,±2,±3}. P. 



SU(3) 



(27) 
denotes 

a projection operator, back onto the SU(3) group and 
the sum is over the four spatial "staples", surrounding 
U^. After extensive studies we employ the parame- 
ter values Nape = 50 for the number of smearing iter- 
ations and the weight factor a = 2.0. For the projec- 
tion operator we somewhat deviate from Ref. [6] where 
we maximized Re tr {A^Pgu^ (A)}, iterating over SU (2) 
subgroups. Instead we define, 



A' = 



VaFa 

P su{3) (A) = A'det(A')- 1/3 . 



(28) 
(29) 



The inverse square root is calculated in the (orthog- 
onal) eigenbasis of A^A, where we take the positive 
root of the respective (positive) eigenvalues. In gen- 
eral there are three possible choices for the phase cor- 
rection, Eq. (29). We take the one that is closest to 



unity. Note that this construction guarantees (except 
in singular cases that we never encountered in numer- 
ical simulations) [Psu(3)(A)]^ 1 = [Psu(3)(A)Y as well 
as gauge covariance: UPsu(3)(A)V = Psu(3)(U AV) for 
U,VeSU(3). 

Subsequently, we construct the spatial transporters Vt 
by calculating products of the APE smeared links along 
paths that stick to the direct connection between quark 
and antiquark as closely as possible. In this way, the 
overlap between creation operator and physical QQ state 
is vastly enhanced. 

We improve the overlap of our mesonic operators 
with the static-light ground state by applying Wuppertal 
smearing [46], 



,(n+l) _ 



1 



±3 



1 + 6(5 



4 n) + *£ u x% 



(30) 



j=±i 



to light quark fields 0, where we set 8 = 4 and replace 
U x ,j by the APE smeared links as detailed above. We 
then employ the linear combination </>( 20 ) — 6.63 23</>( 40 ) + 
7.2604(/>( 5 °) as our smearing function. Note that we are 
calculating local-local all-to-all propagators to which we 
can subsequently apply Wuppertal smearing. 

Best results are obtained by using smeared-local quark 
propagators. For a single static-light meson positivity of 
the coefficients in the spectral decomposition is not guar- 
anteed. Neither do we recover positivity for the bound 
BB system as the source is smeared at position while 
the sink is smeared at position r: the wave function is 
not symmetrized with respect to <-> r. 



2. Static quark action 

One problem in simulations with static sources is the 
rapid exponential decay of the associated Green functions 
with Euclidean time. One of the reasons for this is a large 
static quark self energy contribution which, to leading 
order in perturbation theory, reads, 5m w cCf &l a -1 , 
with a constant c w 1.587956. This contribution obvi- 
ously diverges with a -1 . 

There is however some freedom in the choice of the 
static action, i.e. in the choice of a lattice discretization of 
D4 within Eq. (2), as long as the action remains localized 
and converges towards the continuum action in the limit 
a — > 0. This choice will affect the lattice definition of the 
Schwingcr line Eq. (6). One possible such discretization 
reads, 



DiQ x = a 1 (Q x 



' ^x-ai,iQx-a4, 



(31) 



with 



±3 



U x ,4 = Psu(3) j e£4,4+ £ Ux -i Ux + a i- iU l+ai.j 



(32) 
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where we use e = 0. Note that this procedure is reminis- 
cent of APE smearing, Eq. (27), but with a sum over six 
rather than over four staples. 

The Schwinger line that appears within the corre- 
sponding static propagator can now be written as, 



t/a-l 

U x (t,0)=T Jl F (x>To3>4 . 



(33) 



T denotes the time ordering operator. 

The "extended" temporal links correspond to introduc- 
ing "form factors" in perturbation theory [47] : to leading 
order, replacing U Xi ± by U x> i with the (optimal) weight 
e = is equivalent to multiplying the self energy by a fac- 
tor w 1/2.94: c— 7r/3. The signal is exponentially 
improved in t, while the absolute noise approximately 
maintains its level. Since the self energy cancels from 
energy differences as well as from the sum of 2toq plus 
energy levels, the physics of string breaking remains un- 
affected. Only at small distances, we encounter different 
lattice terms, which (being artefacts of the discretization) 
do not alter the continuum limit. Fat temporal links can 
also influence the ground state overlaps: our impression 
is that they help to improve the situation further. 

One can define a tree level improved lattice dis- 
tance [29, 48], 



[l + 0(a 2 )] 



1 



(r/a) 



(34) 



where [1/R]l — ► 1/R for R — > oo denotes the tree level 
[0(a s )] lattice propagator in an appropriate normaliza- 
tion. It can easily be shown that replacing thin temporal 
links by fat temporal links, Eq. (32), only affects [1/R]l 
at distances R £ {0,£}. Amusingly, iterating the "fat- 
tening" n times will leave the tree-level expressions at 
all distances R\ + R2 + R3 > n invariant. Note that 
= l/(2c) « 0.31a(0.92a) ^ and a w 0.92a(1.37a) ^ a 
for standard (fat) temporal links. In principle, one could 
further fatten temporal links but in view of increasing 
the small-r distortions and of the reduction of the static 
energy already achieved, to one third of its original value, 
we refrain from doing so. Throughout this paper we will 
plot all r-dependent data as functions of r, thus removing 
the short-distance lattice direction-dependence to 0(a s ). 

In the present context, fat temporal links have first 
been employed in Ref. [49]. In this case, more refined 
actions were implemented, utilising all non-intersecting 
paths that can be constructed within an elementary hy- 
percube. Other studies employing similar techniques can 
be found in Refs. [50-53]. 

In Figure 1 we plot effective masses, 



m x , s(t) 



In 



Cx(t) 
C x {t + a)' 



(35) 
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FIG. 1: Comparison of effective masses of static-light corre- 
lation functions, obtained employing static actions with and 
without fat temporal links. The wave function has been opti- 
mized to yield best ground state overlap for the fat link static 
action. SET and HPA have been applied in both cases. 



The curves correspond to one- and two-exponential fits 
to the fat and thin link data, respectively. Note that the 
use of SET and HPA techniques (described in Sec. Ill C 
below) was essential for achieving the high signal qual- 
ity. The new static action shifts the mass by an amount, 
5m = 0.174(7) a" 1 . The absolute statistical errors of the 
correlation function increase somewhat, however, the rel- 
ative statistical errors arc reduced, in particular at large 
times as the signal falls off less steeply. This in turn re- 
sults in a reduction of the error of the effective masses, in 
particular at large t. The Figure also illustrates that we 
are able to achieve an excellent overlap with the physical 
ground state. 

In Figure 2 we compare the static potentials as calcu- 
lated from the Wilson loop operator CQQ{t) alone, with 
and without fat time links. These potentials have been 
obtained from single exponential fits with t m i n /a = 5. 
The curve corresponds to a funnel-type r-dependence 
with fit range 4a < r < 13a. String breaking is ex- 
pected around r — r c rs 15 a. However, this is not seen 
in the data. The error band represents our expectation 
for the thin link potential, obtained from the respective 
static-light mass shift i5m, as determined above. We find 
consistency. 

The extended static action leaves the absolute errors 
of the correlation functions basically constant but still 
improves the signal exponentially. As a result, the effec- 
tive mass errors are reduced impressively, by a factor of 
about five throughout. 7 . 



of static-light mesons with and without fat time links. 
Cs{t) stands for the static- light correlation function. 



7 At first sight, the comparatively modest improvement of the 
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FIG. 2: Comparison of the QQ potential in lattice schemes 
with and without fat temporal links, in both cases for 
tmin/a = 5. String breaking is expected to take place around 
r = r c « 15 a but this is not visible from the Wilson loop data 
alone. The curve represents a funnel fit to the fat link data 
and the error band is this parametrization shifted upwards by 
the amount 28m, where aSm — 0.174(7). 



For completeness we introduce these three techniques 
(TEA, SET and HPA) and their implementation in the 
following subsections. We conclude by comparing nu- 
merical data obtained by use of combinations of these 
methods. 



1 . Truncated eigenmode approach ( TEA ) 

The fcrmionic propagator is the inverse of the Wilson 
Dirac matrix M of Eq. (8). However, M is not Hermitian 
which is why we define, 



Q = J 5 M. 



(36) 



The relation M> = 75M75 implies Hermiticity of Q = 
. We calculate the smallest n = 200 (real) eigenvalues 
qi and corresponding orthonormal eigenvectors |itj), i = 
l,...,n, 



Q\ut 



qi\Ui 



= Si 



(37) 



This is done by means of the parallel implicitly 
restarted Arnoldi method (IRAM) with Chebychev ac- 
celeration [26], using the PARPACK library [55]. We 
can now approximate, 



C. All-to-all propagator techniques 

The requirement of quark propagators from all source 
to all sink locations is obvious: within the Cbb element 
of the correlation matrix, Eq. (23), we encounter light 
quark propagators starting from different source posi- 
tions. Moreover, we can reduce the notorious noise levels 
of disconnected (and of some connected) diagrams; all- 
to-all propagators allow for the full exploitation of trans- 
lational invariance, increasing the accuracy of the entire 
correlation matrix. 

Since the propagator M" 1 , Eq. (8), has 12V x 12V 
components (in our case V = 24 3 x 40), direct evalua- 
tion would be prohibitively expensive, both in terms of 
computer time and of memory. However, the correct re- 
sult can also be obtained by combining the truncated 
eigenmode approach [1, 26, 54] (TEA) with stochastic 
estimator techniques (SET). Where possible, we improve 
both, the convergence of TEA and the statistical errors of 
SET, by employing the hopping parameter acceleration 
(HPA). 



static-light data seems to be in contradiction to the very signifi- 
cant effects observed in Ref. [52], However, a closer inspection re- 
veals that without employing our additional improvement meth- 
ods, i.e. averaging over all lattice points by means of stochastic 
estimates and employing HPA, the gain factor from using the 
extended static action would have been larger: the signal over 
noise improvement appears to saturate, after adding more and 
more tricks. 



E 



(38) 



Obviously, M _1 = Q _1 75- With this truncated eigen- 
mode approach (TEA), we have reduced a 12V x 12V^ 
problem to a 12T^ x 12n problem. 

Another nice feature of this procedure is that there 
is no critical slowing down of the algorithm as the quark 
mass is reduced. However, the difference between the left 
and right hand sides of Eq. (38) is systematic. One can in 
principle estimate this bias from the convergence proper- 
ties under variation of n [26]. In the present context we 
will render the result exact by stochastically estimating 
the remainder, replacing the systematic error by a statis- 
tical uncertainty. For this purpose, it is useful to define 
the projection operator V n , onto the basis spanned by 
the first n eigenvectors, 



Vn = 



n 

E 

i=l 



Ui){Ui 



(39) 



Note that for n smaller than the rank of Q, this basis 
is truncated and hence incomplete: V n 7^ 1. However, 
[Q)V n ] — 0. We can also define the projector onto the 
orthogonal subspace, \ — V n - 



2. Stochastic estimator techniques (SET) 

Stochastic estimator techniques have been applied by 
various groups in the past [2, 27, 56-62]. We introduce 
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the following notation, 



1 N 
N ^ 



(40) 



N denotes the number of "stochastic estimates". Let 
\rf), i = 1, . . . , N be random vectors with the properties, 

(41) 
(42) 



= O(l/VA0 



\r))(rj\ = 1 + 0(1/VN). 



These requirements are for instance met if the 12V com- 
ponents are numbers e*^, with the uncorrelated phases 
4> E {±7r/4, ±37r/4} randomly selected. We employ such 
a complex Z2 noise, where our random vectors take val- 
ues over the entire four- volume, flavour and colour. 
If we solve the linear system, 



Q\s l ) = W), 



(43) 



for \s l ) then, for large N, we can substitute [Eq. (42)], 



\s)(v\- 



(44) 



Note that in our study we actually invert A = M^M = 
Q 2 and then obtain Q^ 1 by multiplying the solution 
with Q. This allows for more flexibility: for instance 
the Roma smearing technique [63] , which amounts to the 
replacement Q^ 1 1— > A^ 1 within hadronic Green func- 
tions, can readily be implemented. It can then be shown 
by means of spectral decompositions that in many cases 
the ground state mass remains unaffected [63]. In the 
present context we have made use of this method, in addi- 
tion to standard smeared-smeared and smeared-local cor- 
relation functions, within the optimization procedure of 
the static-light creation operator. Unfortunately, Roma- 
smearing turns out not to be applicable to the BB <-» QQ 
mixing problem. 

The sparse linear system of Eq. (43) is solved by means 
of the BiCGstab2 algorithm [64]. Unlike in Eq. (38) 
where the bias was systematic, the difference between 
the approximation of Eq. (44) and the exact result is 
purely statistical and reduces like 1/^/N. In order to 
limit the computational effort, N should not be chosen 
overly large. However, the noise level from SET should 
at least match the one from the (finite) sampling of gauge 
configurations. In general, the optimal balance in both 
samplings will also depend on the observable in question 
and on the methods employed. 

We can estimate the difference between the TEA ap- 
proximation and the true result by means of SET. The 
smaller this difference, the smaller the statistical errors 
will be that are introduced by SET. Hence TEA can be 
employed to reduce the variance of SET. We project the 
right hand side of Eq. (44) into the subspace which is 
orthogonal to the TEA eigenvectors: 



(1 - V n )Q-\l - 7>„) « ^ £(1 - V n )W){rf\{l V n ). 

(45) 



In practice this is done by calculating and storing, 

n 

\s j ) = |^>-£k>W>, (46) 

n 

W) = |V)-£k)W>- (47) 

Then, 



i=i 



M~ 



n N 

J2\u i )q7 1 (u i \ + -J2\~s j )(n 



j=i 



75 . (48) 



Note that there is no systematic error on finite-A ap- 
proximants but only a statistical 0(l/\/N) uncertainty. 
In the present context we found N — 50, combined with 
n = 200, to suffice for calculating Cqb^) and static- 
light meson correlators. Within the two diagrams con- 
tributing to Cb b {£) i it is necessary to choose two inde- 
pendent random sources as in either case there exist the 
same two possibilities of connecting sources with sinks. 
In these cases we calculate the SET corrections for the 
two respective propagators independently, with N = 25 
random sources each. Subsequently, we interchange the 
two sets of random sources to increase the statistics at 
little computational overhead. Such "recycling" has been 
pioneered by the Dublin group [65] . 

As long as we are only interested in using SET to re- 
move the bias from TEA for a fixed n we can in principle 
solve Eq. (43) within the orthogonal subspace only, sub- 
stituting the random sources \r} 3 ) on the right hand side 
with |rp). We attempted this but found no advantage 
in terms of real cost in computer time. This of course 
might change at smaller quark masses or with different 
light quark actions. However, having random source so- 
lutions at our disposal that are independent of the TEA 
allows for more flexibility. For instance, not all physical 
states will be dominated by the lowest lying eigenmodes 
of Q = 75M. In particular, the TEA contribution to 
C B °§(t) turned out to be tiny, such that in the end we 
reduced the cost to compute C B °g(t), by employing a 
stand-alone SET. 

Needless to say that, once we have calculated all-to- 
all propagators, Wuppertal smearing, Eq. (30), can be 
implemented. In principle, one could even variationally 
optimize the smearing function [66] , for instance after fix- 
ing to Coulomb gauge. However, our smearing function 
turned out to be already so highly optimized that further 
gain was too hard to achieve. 



3. Hopping parameter acceleration (HP A) of TEA and SET 

The main motivation of complementing SET with TEA 
is to reduce the signal that needs estimation and hence 
the stochastic errors. One might ask if it is possible to 
further facilitate the low eigenvalue dominance, acceler- 
ating the convergence of TEA (and of SET). This is in- 
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deed possible by applying what we call the hopping pa- 
rameter acceleration (HPA). 

We rewrite the fermionic matrix Eq. (8) as, 



M = 1 — kD. 



(49) 



For sufficiently small hopping parameter values k < n c , 
one can expand, 



AT 



oo k— 1 

1 = J2( kD Y = Y.( kD Y + (.KD) k M-\ (50) 

j=0 i=0 



where k > 1. The idea now is that for distances between 
source and sink that are bigger than k lattice spacings 
the first term on the right hand side does not contribute. 
This can readily be seen as follows: as D only connects 
nearest space-time neighbours [and 1 = (kD)° only con- 
tains diagonal entries] , the sum vanishes within elements 
M~y if J2fj, \ x n~Vn\/ a > k > 1. M~y can be replaced by 
[{nD) k M- l ] xy . This means that = [(nDfQ' 1 }^. 
With M\ri) = l-ii\ri), {k\M — Hi{k\, where fii are the 
eigenvalues of M, Eq. (38) can be substituted by, 



Qxy 



M *y^ « E^X 1 - Mi)V 1 ^b)75, (51) 



where again k is smaller or equal to the number of links 
separating source from sink. For large k, contributions 
from big eigenvalues of M are suppressed and the ex- 
pression becomes all the more dominated by low lying 
eigenmodes. Hopefully, the dominance in terms of low 
eigenmodes of M will then also apply to low eigenmodes 
of Q. 

In fact we do not only find HPA to improve TEA but 
the main effect of HPA is with respect to SET 8 : the can- 
cellation of stochastic noise is accelerated if the number of 
contributions to the stochastic average is reduced. Short- 
distance noise is accompanied by larger amplitudes than 
large-distance noise and hence its cancellation requires 
a comparatively larger number of noise vectors. HPA 
explicitly eliminates such short-distance contributions. 

The benefit from HPA will increase with larger tem- 
poral or spatial distances. Unfortunately in the limit of 
light quark masses, as k approaches k c , the quark prop- 
agator will decay less rapidly with the distance and the 
explicit treatment of the first few terms within the hop- 
ping parameter expansion will have less of an effect. This 
is also obvious from the reduced convergence of the hop- 
ping parameter expansion at small quark masses. In this 
case we would however expect a better convergence of 
the TEA contribution in the first place. 



Note that variance reduction methods, that make use of the hop- 
ping parameter expansion, have been pioneered by the Kentucky 
group [60], in a different setting (sec also [27]). 



Note that HPA exploits the ultra-locality of the Wilson 
action and does not generalize for instance to the Neu- 
berger action [67]. Again, this might be compensated for, 
by a faster convergent TEA approximation, due to the 
improved chiral properties of the chiral actions, in par- 
ticular at small quark masses. In contrast, the "dilution" 
method advocated in Ref. [65] will still be applicable in 
a setting with chiral fermions, reducing the variance of 
SET for the very same reasons as HPA does. 

By applying HPA to the whole matrix, cf. Eq. (48), we 
exploit both effects, the improvement of the low eigen- 
value dominance and the variance reduction of SET: 



M~ 



1 N 

i| + a£ 



\~S J )(V J \ 



3 = 1 



75, 



(52) 

where, as above, ka is the lattice-distance between source 
and sink. Again, the above equation is exact up to sta- 
tistical 0(1/ \/N) corrections. 

Unfortunately, due to the size of our smearing function, 
we cannot employ HPA for propagators along spatial sep- 
arations, i.e. within Cqb(£) or within C c ^(t). However, 
Cgg(t) benefits from this technique as do static- light 
correlation functions. One way of extending it to the 
before-mentioned elements is to cut off the radius of the 
smearing function. Such a cut-off, in conjunction with an 
iterative smearing method, is hard to implement. One 
way out would be to work in Coulomb gauge with fixed 
weight-factors [66]. Another possibility is the implemen- 
tation of a rotationally non-invariant smearing function. 
But in this case it turned out to be difficult to sustain an 
acceptable ground state overlap. 



4. Comparative study of SET, TEA and HPA 

We demonstrate the impact of the above methods for 
the example of the static-light meson mass, in the sce- 
nario of the fat link static action described in Sec. Ill B 2. 
We also verify the potential of HPA for the example of 
C/Qsit), however, without smearing (see above). 

In Figure 3 we display effective masses Eq. (35), ob- 
tained with SET alone as well as with HPA SET, and 
after explicitly calculating the contribution from the first 
200 eigenmodes (TEA). Note that the ordinate covers a 
huge t-range, up to a distance t w 1.2 fm. Also note 
the magnified scale of the abscissa, covering the window 
0.9 GeV < mfl :C ff < 1.26 GeV. In particular at large 
times, HPA impressively reduces the errors and TEA re- 
sults in some additional improvement. This is quantified 
in Figure 4, where we display the respective statistical 
errors themselves. Note the logarithmic scale. For in- 
stance, at t = 8a HPA reduces the SET error to about 
one third of its original value while TEA yields another 
w 20 % reduction. Since the effective mass is approx- 
imately independent of t the absolute errors displayed 
are proportional to the relative errors which (as is obvi- 
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FIG. 3: Effective static-light masses, obtained with SET 
alone, with HPA SET and with all three methods combined. 
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FIG. 4: The errors of effective static-light masses, obtained 
with SET alone, with HPA SET and with all three methods 
combined. 



ous from the Figure) grow exponentially with t. Fortu- 
nately, both, HPA and TEA reduce the amplitude and 
the exponent governing this error increase, resulting in 
an exponential improvement at large times. 

In Figure 5 we see that the static-light correlation func- 
tion at large times becomes dominated by TEA. Without 
HPA this dominance seems to be achieved earlier than 
with HPA: the sea quark mass still appears too heavy for 
HPA to significantly enhance the low eigenvalue domi- 
nance of Q = 75M, for the correlation function in ques- 
tion. However, the seemingly perfect agreement for the 
stand-alone TEA case of the displayed ratio at large t 
with one is largely accidental. Increasing the number 



FIG. 5: The relative magnitude of the TEA contribution 
within the static-light correlation function, with and without 
HPA. 



of eigenvectors on one configuration revealed that TEA 
tends to overshoot the exact result, before converging 
towards it. HPA reduces this tendency. Note that HPA 
cannot be applied to the SET part alone, within the SET 
plus TEA combination. 

In the HPA case TEA also results in an impressive re- 
duction of the signal that remains to be estimated. How- 
ever, a comparison between SET plus HPA and SET plus 
HPA plus TEA reveals that the additional error reduc- 
tion due to TEA is only moderate. After HPA the SET 
error is already at the level of the statistical fluctuations 
between gauge configurations and of a comparable size 
to the (non-stochastical) TEA error. In this situation, 
substituting part of one signal by the other leaves the 
resulting statistical error largely unaffected. This would 
have been different for a larger statistical sample or at 
smaller sea quark masses. 

Finally, we wish to investigate the effect of HPA on 
correlators as a function of spatial source-sink separa- 
tions. This is done for the example of Cqb^), without 
smearing 9 . 

In Figure 6 we show the local-local Cqb matrix el- 
ement at fixed t — 4a as a function of r, both for 
SET alone and for HPA SET. All distances are along 
the spatial diagonal, i.e. when increasing F by « \/3a, 
the exponent k increases by three units. At our largest 



9 Our smearing function includes one contribution with 50 Wup- 
pertal smearing iterations which would mean that the exponent 
k has to be smaller than the source-sink link distance minus 49. 
Even for separations along a spatial diagonal, HPA would only 
be applicable for distances much larger than half the lattice ex- 
tent. When using local sources and sinks we do not encounter 
such restrictions. 
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FIG. 6: The effect of HPA on SET for the example of Cqb 
at t = 4a, as a function of the distance r. The normalization 
is arbitrary and we have employed local sources and sinks. 



separation r/a w Ha/3 « 19.05 we have fc = 33! 
The error reduction factors turn out to be fairly time- 
independent. At the largest distance accessible without 
HPA r/a « 8%/3 « 13.9, the error reduction is almost 
five-fold. Potentially, this should be even more impres- 
sive for C"£g , within which two spatial light quark prop- 
agators appear. Unfortunately, the method cannot be 
combined with our present smearing function, which is 
essential for accessing the physical ground states. 



D. Notation and spectral decomposition 

In order to set the stage for the interpretation of our 
numerical data, we detail in this section our notations 
in connection with the spectral decompositions for the 
different matrix elements, Eq. (23). We will assume an 
infinite extent in time direction of the lattice and asymp- 
totic behaviour of all correlators. We often suppress the 
distance dependence, r, from the expressions. 

We define BiBj pair creation operators Bij. These cre- 
ate a light antiquark of flavour j and a quark of flavour 
i, besides the static sources. The states created by Bij, 
i =/= j, constitute a subset of the flavour non-singlet 
("7 = 1") states. The remaining nf — 1 I = 1 states are 
given by traceless linear combinations of the Bij diagonal 
elements. In addition, a flavour singlet (I — 0) creation 
operator, B s = -^=Y^ili^Uj can be constructed. B a 

represents one of the nj — 1 members of the class of I = 1 
operators. 

Flavour singlet states are created both, by B s and by 
Q, the operator that only contains static quark-antiquark 
spinors. We remind the reader of the definitions Eq. (12) 
and Eq. (16), \Q) = Q|0), \B) = S s |0), with vacuum 



state |0). Obviously, (Q\0) = (B\0) = 0. We denote the 
orthonormal eigenstates in the flavour singlet sector as 
|n), n = 1,2, . . ., with energies E n +i > E n . 

For the labelling of the 1=1 sector, we follow the con- 
vention of Eq. (17) and define, \B a ) = £? a |0). All flavour 
non-singlet eigenstates share the same energy spectrum, 
E%. We label eigenstates of the 1=1 sector, within the 
class of states with energy E%> by \n a ), n = 1, 2, 

Just for annotation at this point: while the 1=1 states 
decouple from the 1 = states there will still be mixing 
between 1=1 states and two meson states, containing 
for instance a QQ plus a ir. A calculation of the I = 
1 correlation matrix, analogous to the 1 = sector, is 
beyond the scope of the present paper. 

Note that we smear the B r Bo creation operator at 
position while the sink is smeared at r. This means 
that the creation operator B^ a = Q r ~f5q r <io&j5Qo within 
Eq. (15), where the smearing function $ acts on the 
quark at position 0, is not the Hermitian adjoint of the 
annihilation operator, B s > a ^ ^l\ a - The subscripts "i" 
and "/" stand for initial and final, respectively. Hence, 
strictly speaking, one has to distinguish between {B^\ 
and |-B( a ))- From a practical point of view, the highly 
satisfying ground state dominance of our data renders 
this distinction obsolete. 

We can decompose, 



\Q) = J2(n\Q)\n)=J2Q 

n>l n 

\B i ) = Y / {n\B i )\n)=Y / B: 



(53) 



„,n>, (54) 



n>l 



where 



J2\Qn\' 2 = (Q\Q), (55) 

n 

Y. B CK = (B f \B*), (56) 

n 

Y,Qn B n = J2 B CQn = 0- (57) 



Using these notations, the matrix elements Eq. (23) read 
(neglecting the overall energy off-set 2itiq), 



C QQ (t) 

c QB (t) 
c BB {t) 



?e -E n t 



= ^R C (g;i?;) 



-E n t 



= "£Re(Bl*Bi; 



-E n t 



n 

Note that Re(Q;^) = Re(Q*B^). 



(58) 
(59) 
(60) 
(61) 
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The normalization of our correlation matrix is such 
that C a p(t) > for t — ► oo, i.e. the ground state am- 
plitudes are always positive. However, within Cqb and 
Cbb, excited state amplitudes can be negative. 

For r < r c the ground state |1) will be dominated by a 
|Q)-type component, whereas the first excitation |2) has 
a large \B) contribution 10 . For r > r c this correspon- 
dence will interchange. We view a signal Cqb ^ as an 
"explicit" signature for mixing while a verification of an 
Ei signal in Cqq{1) at r > r c (string decay) or within 
Cbb(0 at r < r c will be referred to as an "implicit" 
mixing effect. 



IV. INVESTIGATION OF INDIVIDUAL 
MATRIX ELEMENTS 

We set the stage for the investigation of the mixing 
problem by comparing individual matrix elements to the- 
oretical expectations. We first discuss Cbb and then 
combinations of different components of the correlation 
matrix, before we attempt to detect implicit mixing ef- 
fects. 



A. Large time asymptotics 

We remark that Cbb (t) contains a disconnected and a 
connected contribution. The disconnected term coincides 
with the 1=1 diagram Eq. (61), and the states it couples 
to are orthogonal to the 1 = sector: any implicit mixing 
can only be mediated through C B °§. This means that at 
< r < r c , 

C BB (t) -> -C B °S(t) (f->oo): (62) 

the connected diagram will dominate at asymptotically 
large t. This is in contrast to the situation at small to 
moderately large times, where the overlaps |2?i| <C I-B2I 
[cf. Eqs. (54) and (60)] warrant a disconnected diagram 
dominance. We shall see below that Eq. (62) in fact turns 
out to be valid for any r > 0, including r > r c . 

We investigate the above expectation in Figure 7, 
where we plot ratios —C B °§(t)/CBB(t) as functions of 
t for a few R = r/a values. Note that Ra = 6a sw ro ~ 
0.5 fm. As expected, at large t, this ratio approaches 
one. Note that r = represents a special case. In this 
limit, the QQ and the BB sectors decouple and Eq. (62) 
does not hold. As a consequence, the asymptotic limit 
is reached faster for r = 2a than for r = a which is ad- 
jacent to r = 0. For r > 2a, the speed of convergence 
decreases again: the gap between the two energy levels 
Ei and E2 reduces as a function of r and hence the limit 



It turns out that (hybrid) excitations of |Q)-nature are energet- 
ically higher than E%. 
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FIG. 7: The ratio —Cbb I Cbb as a function of t for various 
r = Ra. The splines are drawn to guide the eye. 
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FIG. 8: The ratio -C 2 QB I '(CqqCbb) as a function of t for 
various r = Ra. The splines are drawn to guide the eye. 

is approached less and less rapidly in t. At large r the 
signal vanishes in noise, before it can approach unity. 

From Eqs. (58) - (60) it is expected that at asymptot- 
ically large times, for all r > 0, 

\C QB (t)\ 2 -> C QQ (t)C B B(t) (t-00). (63) 

In Figure 8 we verify this expectation for some selected 
distances. The disconnected contribution to Cbb has to 
decay faster than the connected contribution. Otherwise 
the corresponding n t , colour and Dirac factors would be 
incompatible with Eq. (63). Hence Eq. (62) above is not 
only valid at r < r c but for any distance r > 0. 
In particular this means that at large t, 

C% B (t) - -C 00 (t)CgS(t), (64) 
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but is only contained within Cqq. At r > r c the ground state 

FIG. 9: The ratio Cbb /C b 

as a function of t for small r — Ra. ^ s on i v visible in Cbb 



or, diagrammatically, 
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X 








X 





(65) 



again for any r > 0. This limiting behaviour is ap- 
proached faster in time than Eq. (63) above. 

Both sides of Eq. (63) are dominated by the ground 
state contribution |1), which results in an exponential 
decay oc e~ 2Elt . At r < r c , Cqq couples more strongly 
to this term than Cbb- This interchanges at r > r c , 
where the ground state is dominantly contained within 
the Cbb sector. The decoupling of Cqb from the 1=1 
sector implies, E\ < Ef: the 1 = ground state energy 
cannot be larger than the lowest 1=1 energy level, at 
any distance r > 0. 

Finally, we compare C^g(t) to the static-light corre- 
lation function C_b(<). If the i?-mesons at positions 
and r did not interact with each other then the ratio 
C B %(t)/C B (t) would be unity. The I = 1 BB state 
would merely act like the sum of two isolated B mesons. 
We investigate this ratio in Figure 9 and find this sce- 
nario to be valid within our statistical resolution for 
r > 2^/2a 0.23 fm. The increase of this ratio at large 
t for small r < 2\[2a can be attributed to an increased 
overlap of the creation operator with the QQir 1 = 1 
ground state, see also Sec. VD. 



B. Implicit detection of mixing effects 

We investigate whether the actual ground state energy 
level is visible in Cbb^) at r < r c or in CQQ(t) at r > 
r c . For this purpose, we study the large t behaviour of 
effective masses, Eq. (35). The qualitative expectation 
in the absence of any mixing is sketched in Figure 10: 
at r < r c the ground state can only be detected within 
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FIG. 11: Effective masses. The solid symbols are from 
CsB(t), the open symbols from the Wilson loop CQQ(t). The 
horizontal line corresponds to twice the mass of a static-light 
meson. 



Cqq, while at r > r c the ground state energy is given 
by the large t behaviour of Cbb and will not be visible 
from Cqq. Implicit mixing means that the Cqq and 
Cbb effective masses share the same ground state. At 
r < r c the Cqq effective mass is expected to plateau at 
smaller t values than the Cbb effective mass. At r > r c 
the ground state then will become dominated by the \B) 
contribution and hardly be visible in Cqq. 

We display the situation for small r in Figure 11. The 
open symbols, calculated from Wilson loops Cqq, ex- 
hibit good and early plateaus. The solid symbols, which 
correspond to the matrix element Cbb, start out from 
values, similar to the mass of two static-light mesons, 
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FIG. 12: Effective masses from the Wilson loop Cqq. The 
lowest horizontal line is 2ms- The other two horizontal lines 
denote the central values of the respective excited state expec- 
tations from a global linear plus Coulomb fit to the potential 
from 2 < R < 13 Wilson loop data. The open pentagons, 
from Cgg, demonstrate the plateau quality within the Cbb 
sector. 



but then decay towards the respective lower lying states, 
clearly signalling implicit mixing effects! Note that as 
C BB alone only projects onto the 1=1 sector, this effect 
is entirely due to the C C BB contribution, see Figure 7. 

We are also tempted to verify implicit string breaking 
at large r. To this end, in Figure 12 we examine effec- 
tive Wilson loop masses for two distances r > r c as 15a: 
r/a = 10\/3 w 17.32 and r/a = ll\/3 w 19.05. The 
upper two horizontal lines denote the respective plateau 
value expectations from a global linear plus Coulomb fit 
to the potential as obtained from Wilson loops, in the 
region without string breaking, r < 11a < r c . We wish 
to see the data to deviate from this expectation, towards 
the lowest line, that corresponds to 2ms. We conclude 
that we see no implicit indications of string breaking in 
the Wilson loop data. 

We also demonstrate the quality of the effective mass 
plateau within the Cbb sector in Figure 12 (open pen- 
tagons). Here, in the interest of small error bars, we ne- 
glected the C B ° B correction to the masses. As expected, 
the C BB effective mass data agree with twice the static- 
light mass. We remark that we also find E\ as Ef within 
errors, in agreement with the level ordering expectation, 
Ei < Ef. 

The absence of an indication of implicit string break- 
ing at large r is no surprise since in a study of adjoint 
potentials in 2 + 1 dimensional SU(2) gauge theory [22], 
this was only seen at physical times much bigger than 
ours. We shall address the question, at what t values im- 
plicit string breaking should become visible, in Sec. V C 
below. We will see in Sec. VA, where we also study ex- 



plicit mixing effects, that mixing is indeed much smaller 
at r > r r than at r < r r . 



V. RESULTS 

We present our analysis and results on the 1 = mix- 
ing angle and energy levels. We then discuss string break- 
ing as well as transition rates, before we address the short 
distance behaviour of the energy levels, both within the 
1 = and the 1=1 sectors. As we only work at a fixed 
value of the lattice spacing, it is convenient to display all 
results in this section in lattice units, a sa 0.083 fm, i.e. 
a" 1 w 2.37 GeV. 



A. The mixing analysis 

Our creation and annihilation operators are highly op- 
timized, such that the overlaps Q n = (n\Q) and B n = 
(n\B) are close to zero for n > 3. Hence, we base our 
analysis on the simplified mixing scenario, 



\Q)=a Q (cos0|l>-sin0|2)), 
\B) = a B (sin0|l) +cos0|2)), 



(66) 
(67) 



truncating Eqs. (53) and (54) after n = 2. In what fol- 
lows, we abbreviate, c$ = cos0 and sq = sin0. The 
identification, Q\ = aQCg, Q 2 = —Qqsq, B\ = a B sg and 
B2 = CLsCe above guarantees that (Q\B) = 0, as well 
as positivity of all correlation matrix elements for large 
Euclidean times. 

The ansatz Eqs. (66) - (67) implies that, 

C QQ {t) = a 2 Q [c 2 e exp(-E 1 t)+s 2 e exp(-E 2 t)], (68) 
C BB (t) = a 2 B [s 2 e e^(-E 1 t) + c 2 8 cxp{-E 2 t)], (69) 
Cqb(*) = a Q a B sece [exp {—E\t) - exp (~E 2 t)] .(70) 

We also determine the correlation functions at t = 
0. This enables us to implement the normalization, 
Cqq(0) = Cbb{0) = 1. In this case, the ciq and as val- 
ues, fitted at large t, can be interpreted as the overlaps of 
our respective trial wave functions with the n < 2 eigen- 
state sector, with optimal value, clq = a 2 B = 1. While 

o-q = 1, this is n °t necessarily so for a 2 B as here addi- 
tional exponentials can come in with negative weight. 

We attempted to model corrections to Eqs. (68) - (70), 
by adding additional exponentials to our fits. The over- 
laps (n\B) for n > 3 were so tiny that, except at r < 2a, 
we were unable to detect any additional masses in the 
Cbb channel. It was however possible to add an addi- 
tional excitation to the \Q) channel. This exponential 
then also couples to the Cqb element. Results of such 
eight parameter fits were very compatible with those of 
the simultaneous five parameter fits introduced above, 
with parameters 0, aQ, a#, E\ and E 2 - However, these 
eight parameter fits were not stable at all distances. In 
contrast, the five parameter fits turned out to be very 
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robust, such that the results presented here are based on 
the parametrization Eqs. (68)- (70). We also attempted 
six parameter fits, allowing 9 to take two different values 
within Eq. (66) and Eq. (67). This however did not im- 
prove the x 2 /Ndf qualities and the two 0-angles turned 
out to agree within errors. We conclude that the mixing 
scenario Eqs. (68)-(70) is preferred by the data. 



^QQ.min — 4a, t QB ,min 
^QQ, min = tQB ,min 
,min 

tQQ, min = 6a, tQB, min 
t<QQ, min = 6a, tQB, min 

We determined the correlation matrix elements C a f} it) 
for < t < 10. However, the quality of the Cqq and 
Cqb data only allowed us to use t max = 9 for r > 12. 
For the discussion of the individual matrix elements pre- 
sented in the previous sections, we calculated jackknifc 
and bootstrap errors. Both showed consistent results. 
Hence, in this more complicated mixing analysis we re- 
strict ourselves to the jackknife method. 

Prior to the fits we transformed the data: 

C af s(t)^C a0 (t)/C 2 B (t). (76) 

This automatically normalizes the energy levels with re- 
spect to 2m b, removing the self energy 2toq. In prin- 
ciple, this procedure bears the risk of introducing addi- 
tional energy levels that are present in the static-light 
correlation function but not within the QQ/BB system. 
Figure 3 confirms nicely that at t > ft no such levels are 
statistically detectable. Cb(0), however, turns out to be 
about 2 % larger than an exponential extrapolation down 
from t > a suggests. We correct for this in the calcula- 
tion of the overlaps oq and as- The results for all fit 
parameters are summarized in Table I. 

9 is the mixing angle of the physical eigenstates |1), |2) 
with respect to the Fock basis \Q),\B), used in the simu- 
lation. We invert Eqs. (66) - (67), adapting the normal- 
ization \Q) i ► a^\Q), \B) i ► a^l-B): 

|1) = cos9\Q) + sin0|.B), (77) 
|2) = -sin0|Q) + cos0|B). (78) 

In our tap > t Q /3,min analysis we effectively encounter 
this idealized picture and truncate the eigenbasis at n — 
2, which is supported by the data. However, we also 
truncate the Fock basis after the QqqQ sector. In general, 
the physical eigenstates will receive contributions from 
higher Fock states too and hence there will be a (slight) 
model dependence in the determination of 0. 



For each r we carefully checked the quality of the fits 
and the stability of the parameter values with respect to 
variations of the fit range. For each of the three matrix 
elements we determined a i Q/ 3, m i n value, as turned out to 
be much closer to unity than ciq and we found, iQQ, m i n > 
tQB, min > tBB,min- The fit ranges that we employed in 
our final analysis are, 



= 2a for r > 15.5a, (71) 

= 2a for 15.5a > r > 13a, (72) 

= 2a for 13a > r > 12a, (73) 

= 2a for 12a > r > 8a, (74) 

= 3a for r < 8a. (75) 

I 

The limit r = represents a special case. In this limit, 



cxO- (go7(l -74)<Zo)c/ = 0, (79) 

since (q r yiq)u = (qjij4q)u = 0: the two 1 = eigen- 
states decouple and hence 9 = 0. Moreover, {W(0, t))u — 
const, which means that the vacuum is the ground state. 
Q and Q annihilate: £a(0) — 2ms = -2mg. The 
ground state overlap at this point is ag = 1, by def- 
inition, and Cqb is undefined. At r = we perform 
a two exponential fit for t > a to CsB{t), with expo- 
nents Eg and E^. The lower mass (without subtract- 
ing Ims) is Eg = 0.32(17)a _1 but with tiny overlap: 
a 2 g = 0.0045(15). This mass should coincide with the 
mass of two interacting pions or with the scalar /o mass, 
see also Sec. VD below. We have 2m w = 0.55(l)a _1 for 
two non-interacting pions while the mass of the lightest 
scalar is [6] m ++ = 0.71(6)a _1 . Both values are com- 
patible within two standard deviations with the fitted E g 
value above. 

B. String breaking 

In Figure 13 we plot the two energy levels, normal- 
ized with respect to 2ms, as a function of the 0{a s ) 
improved lattice distance r, Eq. (34). Note that string 
breaking takes place at a distance r c « 15a ~ 1.25 fm. 
The implications with respect to the nj = 2+1 QCD situ- 
ation with realistic quark masses are discussed in Sec. VI 
below. The curve corresponds to the three-parameter fit, 

E 1 (r)=V + ar-e/r, (80) 

with fit range 0.2 fm < r < 0.9 fm. For the normalization 
we find, V = 2m B — 0.509(8)a~ 1 while string tension and 
Coulomb coefficient are respectively, 

= 0.1888(29) a" 1 , (81) 
e = 0.362(16). (82) 



I 

= 3a, tsB,min 
= 4a, t B B,min 
= 4a, t B B,min 
= 5a, ts_B,min 
= 5a, tBB,m\n 
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TABLE I: The energy levels Ei and E2, as well as the mixing angle 6, the overlaps ciq and as and the transition rate g, defined 
in Eq. (94). Note that at r = 0, E\ = and the ground state energy is, a(E g — 2ms) = —0.59(17). E2 denotes the first 
excitation above this state. 
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FIG. 13: The two energy levels, as a function of r, normal- 
ized with respect to 2ms (horizontal line). The curve cor- 
responds to the three parameter fit to Ei(r), Eqs. (80)-(82), 
for 0.2 fm < r < 0.9 fm < r c . 
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FIG. 14: The same as Figure 13, for the string breaking re- 
gion. 



The fit implies a Summer parameter 
2 dE^r) 



'o 



dr 



1.65, 



of 



r = 6.009(53)a«0.5fm, 



(83) 



(84) 



which we use to translate the lattice scale a into physical 
units. 

On the scale of Figure 13, the energy gap AE C = 
min r [_E 2 (f*) — E\(r)\ is barely visible. Therefore, we en- 
large the string breaking region in Figure 14. We define 



the string breaking distance as the distance where the 
energy gap is minimal: i?2( r c) — Ei(r c ) — AE C . 
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FIG. 15: The mixing angle 9, as a function of r. The curve 
corresponds to the parametrization Eqs. (85)-(88). 

Not only the two energy levels play a role in the mixing 
dynamics but also the mixing angle 9 of Eqs. (77) and 
(78). In Figure 15 we depict 9 as a function of r. For r < 
r c , the overlap Q\ will be larger than B\ and hence 6 < 
7r/4. For r — > oo the QQ content of the ground state will 
vanish and 9 — > ir/2. The Figure reveals that while this 
large r limit is rapidly approached for r > r c , the ground 
state at small r contains a significant BB admixture: 
for instance, sin 2 [6>(8a)] « 0.03. Furthermore, there is a 
"bump" at small r in 9(r) as well as in ^(r), before 9 is 
forced to approach zero at r -> U , where C QB (t) = 0. 
This bump is likely to be related to light meson exchange, 
where in our study m~ x « 4a. 

The curve corresponds to a phenomenological three pa- 
rameter fit to the 0.9 fm ss 11a < r < 19a w 1.6 fm data: 



9(r) - | {arctan [d(r - r a )] - |} + ^ , (85) 
with parameter values, 



(86) 
(87) 



r s = 14.95(12)a, 
d = 2.31(21)a" 1 , 
c = 0.914(6). 



The increase of 6 with respect to r for r ~ r s is given 
by, d9(r)/dr\ r=rs = cd/2 = 0.34(3)7ra _1 . Our distance- 
resolution clearly allows us to resolve the mixing dynam- 
ics at r ~ r c . We enlarge this region in Figure 16. 

Finally, in Figure 17, we investigate the difference 
AE(r) — E2(r) — E\(r) in the string breaking region. 



Note that ss 0.92a. 
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FIG. 16: The same as Figure 15, for the string breaking re- 
gion. 




FIG. 17: The energy gap AE — E- 2 — Ei with (circles) and 
without (squares) mixing. 



The circles represent the results from our mixing analy- 
sis while the squares are extracted from fits to the Wilson 
loops Cqq and the 1 = 1 BB operator Cg s B alone. This 
resembles the situation in the quenched approximation 
where no string breaking or mixing occurs. We perform 
a quadratic fit in the region 14a < r < 16a, 



AE{r) = AE C + b 3 (r - r c f 

The resulting parameter values are, 

r c = 15.00(8) a, 
AE C = 0.0217(9) a" 1 , 
b = 0.325(14) a" 1 . 



(89) 



(90) 
(91) 
(92) 



The position of the minimal energy gap r c = 15.00(8)a 
is in perfect agreement with the value r s = 14.95(12)a 
of Eq. (86), at which 8 = ir/4. Translated into phys- 
ical units we obtain a minimal energy gap, AE C « 
51(3) MeV, and a string breaking distance, 



2.496(26) r w 1.248(13) fm. 



(93) 



The errors quoted are purely statistical and do not con- 
tain the 5 % uncertainty of ro ~ 0.5 fm or the deviation 
of rif = 2 and m < m s from the real QCD situation. 
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FIG. 18: The transition rate g between \B) and \Q) states, as 
a function of f. 



C. Transition rates 

We assume that the elements of our mixing matrix 
only couple to the lowest two QCD cigcnstates within the 
appropriate static-static sector. In this limit, for each r, 
we encounter a quantum mechanical two-state system. 
Our two test wave functions are not QCD eigenstates 
and, therefore, the off-diagonal matrix elements Cqb(£) 
assume non-trivial values. The transition rate, governing 
string fission at r > r c and fusion at r < r c , is given by, 



dC QB (t) 



dt 



;=o V c bb(0)C qq (0) 



(94) 



While in Euclidean time all Fock states eventually de- 
cay into the ground state |1), in Minkowski space-time, 
starting from such a non-eigenstate, results in oscillations 
between the QQ and BB sectors. 

Obviously, our states \Q) and \B) are somewhat pol- 
luted by n > 3 excitations as evidenced by oq =/= 1 and 
as ^ 1. So we have to "wait" for some initial relaxation 
time t m ; n to pass until this equation becomes applica- 
ble. We can easily extract g from our five parameter fits, 
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Eqs. (68) - (70), setting C QQ (0) = a 2 Q and C BB (0) 



g(r) = AE(r) 



sin[20(r) 



(95) 



This quantity is plotted in Figure 18 and the re- 
sulting values are displayed in the last column of Ta- 
ble I. At the string breaking point, 9 w n/4 and hence 
g(r c ) « AE c /2 w 25(2) MeV. This means that g(r) has a 
maximal value of around 320 MeV, at a distance of about 
0.2 fm. The curve is a polynomial and drawn to guide 
the eye. g~ l can be interpreted as the characteristic Eu- 
clidean time scale, governing the decay of one state into 
the other. For instance, at r — 2a we find the (maximal) 
value g~ x w 7.3a and indeed, in Figure 11, at t/a > 7 the 
Cbb effective masses (solid circles) agree with the Cqq 
level (open circles). We also find that the implicit indi- 
cations of mixing effects are most pronounced at exactly 
the distance at which g is largest. 



For small r, g decreases as it has to reach zero at 
r = 0. At large r, 9 approaches it/2 quite rapidly, 
resulting in small g values too: for r > r c we find 
g < AE c /2 ps 0.011 a -1 . This means that detecting the 
(dominantly \B)) ground state of the system from Wil- 
son loop signals alone necessitates distances t = O(100 a). 
Possibly, depending on the statistical accuracy, t = 50 a 
might be sufficient to verify the decay of the Wilson loop 
signal towards the ground state energy. In view of this, 
it is no surprise that in Figure 12 we have been unable to 
verify such implicit string breaking at r > r c from t < 9a 
data. 

It is possible to calculate g directly from the data, 
without any fits. This will be a valuable consistency 
check. For this purpose, the time derivative has to be 
eliminated from Eq. (94) . It is straight forward to derive 
the approximate expression (for a similar result, see e.g. 
Michael [68]), 



g(t) » Z 



C Q B(t)y/C QQ (0)CBB(0) 



Et I a — 1 



C QQ (ja)CBB(t - ja) + i [C QQ {Q)C BB {t) + C QQ (t)C BB (0)] 



Z = 



aAE 



2 tanh(aA£y 2) ' 



(96) 



The Z term originates from replacing a time integral by 
a discrete lattice sum. The required level difference AE 
can be approximated by an effective mass, however, for 
our proof-of-principle calculation we use the AE values 
of Table I, extracted from our five parameter mixing fits. 
The largest such correction for the examples, displayed in 
Figure 19, amounts to a 3.7 % upward shift of the r = 2a 
data. 

If the energy gap AE is large then, within the denom- 
inator, the propagation of the lighter state is strongly 
preferred over that of the heavier state and the transi- 
tion between the two states will take place near the end 
points j w or j m t/a. In this case, unless E s ^> E 2 , 
there will be higher state contaminations and no accurate 
result can be expected. If AE is large then g = sgcgAE 
can also be large. In the derivation of Eq. (96) implicit 
mixing effects are neglected and due to this, at large t, 
there will be corrections, g = g(t) + O(gt). If g is suf- 
ficiently small, then there is a chance of identifying a 
plateau in g(t) from large enough (but not too large) t 
values. 

In Figure 19 we compare g(t) approximants, obtained 
by use of the modified Michael ratio method Eq. (96), 
with our fitted g values (horizontal error bands). We 
find good plateaus and perfect agreement with the fitted 
g values, except at distances r < 6a where implicit mix- 
ing is significant and the linear t behaviour already sets 
in, before the excited state contributions have died out. 
In principle, one could attempt to subtract such linear 
terms. 
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FIG. 19: The modified Michael ratio Eq. (96) for different 
r = Ra values, as a function of t. The horizontal error bands 
correspond to the fitted transition rates. 



The ratio method offers a nice check of consistency. 
Other than this, we see little advantage in calculating g(t) 
over extracting g from a global fit (in t) of the correlation 
matrix elements. In the case of very noisy data such 
fits may turn out impossible but in this case any g(t) 
estimate will be unreliable anyhow. Note that within our 
fit ranges Eqs. (71)-(75), g(t) does not show any sizeable 
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excited state contaminations. 



D. Short distance BB forces 

We now focus on the short distance behaviour, both 
for / = and for / = 1 BB systems, and display the 
respective two lowest lying energy levels in Figure 20. 
Note that the difference between the 1=1 and the 1 = 
ground states is given by our (unrealistically heavy) tt 
mass, m„ ss 640 MeV. The lowest state in the 1 = sec- 
tor corresponds to the conventional QQ potential. We 
already noted the bump in the excited state level, with 
a maximum of about 85(10) MeV, relative to 2ms, at 
a distance of 0.2 - 0.3 fm. Such a bump is not present 
within the 1 = 1 sector, to which the exchange diagram 
Cgg does not contribute. We assume this energy barrier 
to be related to meson exchange. Note that at the dis- 
tance of the maximum, AE(r) ss Im^. Unfortunately, 
in our study we restricted ourselves to one quark mass 
and hence we are unable to investigate the quark mass 
dependence of the height and of the position of this fea- 
ture. 

Within the 1 = 1 situation, we also encounter two 
sectors, namely a QQtt state (with QQ in the £~ repre- 
sentation that is mass-degenerate with E+) and the BB 
state that we label as \B a ). After diagonalization of this 
mixing problem one should be able to identify a mixing 
angle 4> and the two energy levels 12 Ef and E% , in anal- 
ogy to the 1 = system. Again, for r = the off-diagonal 
elements of the corresponding correlation matrix vanish 
and the two sectors decouple. The lower lying state will 
be a single tt, with Q and Q annihilating. In contrast, 
in this limit, Cgg will couple to scalar states. In the 
colour singlet sector this will be a scalar ao meson as 
well as 7T7T scattering states. There will be excitations 
above these mesonic states, corresponding to two light 
quarks, bound to an adjoint static colour source, in anal- 
ogy to pure gauge hybrid potentials where Q and Q do 
not annihilate at r = (gluelumps [34]). Note that in 
the limit r = 0, the 1 = correlation function CBs(i) 
will also couple to both, BB states with qq in a colour 
octet (which we shall call qq lumps) and to colour singlet 
fo/mr states. The latter sector is lighter. 

As we only have the creation operator of the \B a ) state 
B a at our disposal but do not separately investigate the 
QQtt sector (and the mixing between the two sectors), 
we assume the lowest lying 1=1 state (open squares) 
to consist of the ground state potential plus the mass of 
the tt. There will be a (small) correction to this assump- 
tion, due to the interaction energy on a finite lattice. For 
the r > 2^/2a data we are unable to detect this state 



12 Note that at very small distances and/or small sea quark masses 
there will be additional multi- mesonic states between the |l a ) 
ground state and the \2 a ) excitation. 
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FIG. 20: The first two energy levels within the 1 = and 
1 = 1 sectors, at short distances. 
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FIG. 21: The 1 = and 1 = 1 excited state energies, relative 
to 2m b, at short distances, r < ro sa 6a. 



within the Cgg signal, see also Figure 9. For r < 4a we 
apply two exponential fits with i m j n = 2a. The overlap 
with the QQtt ground state turns out to be tiny in these 
three parameter fits, with mixing angles ranging from 
sin 2 (c/>) = 0.0006(22) at F w 3.51a to sin 2 (c/>) = 0.0034(4) 
at r ~ 1.37a. However, the fits are consistent with the 
ground state mass assumption, E\(r) + m n . 

In Figure 21 we focus on the two E2 levels at small r. 
As noted before, there is repulsion in the 1 = sector for 
distances r > 0.25 fm, with a peak value of the energy 
barrier of about 85 MeV. However, at very short range, 
attraction sets in. This has to be so since in the r — * 
limit the first excitation above the vacuum is a flavour 
singlet light quark state, with Q and Q annihilating each 
other. Hence the form at asymptotically short distances 
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will be governed by the perturbative colour singlet po- 
tential. 

Contrary to Ref. [61], we also observe (weak) repulsion 
in the 1 = 1 sector. This difference might be due to a 
bigger overlap of the B a operator used in this previous 
study with the QQtt ground state. However, with our 
operator and statistical accuracy we are able to clearly 
separate the (tiny) QQn pollution from E 2 . From a four 
parameter two exponential fit to the 1=1 operator C BB 
at r — we find, aE\ = 0.394(26), very similar to the 
corresponding 1 = value, aE g = 0.32(17). This might 
indeed be a scalar a meson. The coupling between our 
operator and this state is a\ = 0.010(2). The first exci- 
tation (with which our operator has 99 % overlap) that 
we are able to resolve is E 2 — 2ms = 0.0079(16)a~ 1 (left 
most data point: w 0.92 a). We interpret this as the 
lowest lying lump of a qq state, bound to an adjoint static 
colour source. In this case the short distance interac- 
tion can be identified with the octet potential [34]. As 
argued above, there should be further scattering states 
inbetween the ao level and the qq lump, however, our op- 
erator basis appears to have (almost) vanishing overlap 
with them. 



VI. PHENOMENOLOGICAL IMPLICATIONS 

We discuss a possible extrapolation of our string break- 
ing results to the rif = 2+1 case with realistic light quark 
masses. We also comment on the relevance of the results 
with respect to quarkonium spectroscopy. 



A. Extrapolation to real QCD 

We expect the string breaking distance to decrease 
with the sea quark mass. From the experimental differ- 
ence m Bs — itib = 90(3) MeV, we obtain r c w 2.16 r w 
1.08 fm if we assume invariance of the shape of the QQ 
potential under variation of the sea quark mass. This 
assumption however is rather arbitrary and we wish to 
refine this first very rough estimate. 

A more controlled way is to extrapolate previous re- 
sults of the static-light meson mass [38, 69] and of the 
QQ energy [6] quadratically in the 7r mass. The latter 
extrapolation has already been performed in Ref. [6] . For 
the static- light mass we obtain an upward shift, Am B = 
0.21(10) rg 1 , when replacing our simulated quark mass 
by the physical light quark mass. This direction of change 
is possible since the self energy of the static propagator 
increases with decreasing sea quark mass [70]. The po- 
tential at r c also moves upwards, unsurprisingly by an 
amount that is larger than 2Atob. In combining the 
two extrapolations we obtain r c = (2.27 ± 0.20) r ~ 
(1.13 ± 0.10) fm for rif = 2 light sea quarks, in good 
agreement with the rough phenomenological estimate 
presented above. 



The previous lattice results [6, 38, 69] were obtained 
with a static action that differs from the present one 
where we employ fat temporal links (see Sec. IIIB2), 
however, this change will not affect the string breaking 
distance since, when introducing the fat link action, both 
energy levels are always shifted downwards by the same 
amount. 

We discuss the effect of a third, heavier, sea quark 
flavour. In this case there will be two separate thresh- 
olds, one for the decay into what we call B and B mesons 
and one into B s and B s . It is not a priori clear what 
effect the inclusion of such a third sea quark has on 
the r c position at which the decay into BB sets in. A 
comparison between the rif = and the n/ = 2 sit- 
uations might give some indication. Interpolating the 
rif = static-light masses of Refs. [71, 72] to our 
quark mass, m^/my = 0.704(5), we obtain the value 
m B = 0.540(10) a- 1 at = 6.2 where r = 7.30(4) a. To- 
gether with the potential from Ref. [30] , this corresponds 
to r c = 2.53(8) ro, very consistent with our rif = 2 result, 
Eq. (93), r c = 2.50(3) ro. So we would expect the value, 

r c = (2.27 ± 0.20) r w (1.13 ± 0.10) fm, (97) 

to remain largely unaffected by the addition of the 
strange quark. Note that there are additional systematic 
errors of about 5 % on the scale ro and that we have not 
attempted a continuum limit extrapolation. We expect 
large distance physics like the string breaking scenario 
to remain largely unaffected by charm quark dynamics 
which, however, might influence short distance interac- 
tions. 

In Figure 22 we display our rif = 2, to < m s energy 
levels in physical units. The plotted parametrizations 
are, 



where 



Ei(r) 
Mr) 



9i{r) 
V(r) 



2m B +g 1 {r)V{r)+c 1 , (98) 
2m B + [l-g 2 (r)]V(r) + c 2 , (99) 



^ - - arctan [dj(r - r c )] , (100) 

i 7T 

-e(-- -) + <j(r-r c ). (101) 
V r r r J 



We use the arctan function in the definition of the 
smeared out step functions gi(r), rather than e.g. tanh, 
to allow for a direct comparison with the dependence of 
the mixing angle on r, Eq. (85). Also note that the 
above parametrizations represent only effective descrip- 
tions of the data, within a certain window of distances 
r < 1.6 fm. For instance, E\ does not have the correct 
large distance limit 2m B . The parametrization of E\(r) 
is valid for r > 0.2 fm, while that of E 2 (r) applies to 
r > 0.75 fm. In this latter channel, we encounter a repul- 
sive potential barrier at smaller distances, see Figures 13, 
20 and 21. 

We use the same e w 0.36, ^fa k, 447 MeV as in 
Eqs. (81) - (82) and the r c w 1.25 fm of Eq. (93). We 
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FIG. 22: The energy levels in physical units for nj — 2 at a 
quark mass of slightly less than the strange quark mass (data 
points). The bands represent our rif = 2 + 1 speculation. 



obtain a « -21 MeV and c 2 ~ 31 MeV. Note that 
AE C = C2 — c\ pa 51 MeV. The parameters cfj read, 
di w 2.7 GeV (0.073 fm)" 1 , and d 2 « 4.0 GeV ps 
(0.049 fm) -1 . Note that we parameterized the mixing 
angle in the string breaking region in a similar fashion, 
Eq. (85), with <2ps 5.5 GeV. 

We speculate about the real QCD situation in Fig- 
ure 22 (bands). Besides the above discussed reduction 
of the string breaking distance, we would expect the 
shape of the energy gap AE(r) to depend on the quark 
mass as well. A lighter mass will result in a larger gap 
AE C and a broadened mixing region. We plot the cor- 
responding curves with the arbitrary correction factors, 
Cj h- > 1.5 Cj, di i— > di/1.5, also taking into account an in- 
crease of e pa 0.4 and a reduction of pa 440 MeV [28]. 
There will be a second level crossing around 1m^ s , which 
we also sketch in the Figure. 

Another uncertainty is related to the short distance 
dynamics that we observe. We found an 85 MeV high 
potential barrier within E^ir), at a distance of about 
0.2 fm. Most likely this is related to ir exchange [73, 74]. 
In this case the dimensions, both of the height of the 
barrier and of its position should be provided by the 7r 
mass m^. On one hand, reducing by a factor of four 
down to its physical value could easily move this region 
close to the string breaking distance. On the other hand, 
we would then expect the associated correction to the 
E 2 (r) level around r c to be smaller than 20 MeV, while 
AE C > 50 MeV. 

Finally, we remark that at smaller quark masses ad- 
ditional scattering states will occur between the Ei(r) 
and the E 2 (r) energy levels, at short distances where 
AE(r) > 2m^. It is clear that replacing our qualitative 
Uf = 2 + 1 picture by a truly quantitative understanding 
would require simulations at additional light quark mass 
parameters. 



B. Relation to quarkonium physics 



String breaking provides an intuitive example of a 
strong decay. In addition, static potentials can readily 
be related to quarkonium physics: this can be achieved 
by introducing a phenomenological Born-Oppcnhcimcr 
approximation [75] or, more systematically, within the 
framework of potential NRQCD (pNRQCD) [76]. As 
long as the quarkonium state in question is much lighter 
than the respective strong decay threshold into a pair 
of heavy-light mesons, the ground state energy level 
Ei(r) should, to leading order in the relative QQ veloc- 
ity v <C 1, accurately encapsulate the relevant physics. If 
this is not the case anymore, then additional terms have 
to be added to the pNRQCD Lagrangian, to incorporate 
the QqqQ sector and transitions between the two sec- 
tors: when it comes to strong decays like Y(45) — + BB, 
E 2 {r) and the transition rate/coupling g(r) [or, equiva- 
lently, the mixing angle 9{r)\ are required to describe the 
system, in addition to E\(r). 

For states that are stable against such decay, but which 
are in the vicinity of a threshold, mixing effects will result 
in mass shifts [75, 77]. This has been discussed in some 
detail for instance in Refs. [74, 78, 79] for the recently 
discovered X(3872) charmonium state [80, 81]. In this 
context, our results could hint towards the nature of the 
underlying interaction Hamiltonian and of the qq pair 
creation mechanism that is at work. 

We have demonstrated in Sec. VC that mixing is suf- 
ficiently weak to allow a basis of QQ and QqqQ quark 
model Fock states to be a valid starting point in any such 
analysis [82]. However, not only in the string breaking 
region but also at short distances the ratio g(r)/ AE(r) 
can be sizeable. We find g(r) to be of an O(100 MeV) 
magnitude which is typical for strong decay dynamics. 
g(r) sets out from zero at the origin, increases to about 
320 MeV around r pa 0.2 fm and reduces to AE c /2 pa 25 
MeV, in the string breaking region. The maximum value 
is due to meson exchange and its position coincides with, 
AE(r) pa 2m T pa Ag{r). We would expect this medium 
range g value to somewhat decrease with lighter quark 
masses and g(r c ) to increase. 

In QCD with sea quarks there are not only BB excita- 
tions present but also QQ-gluon hybrid potentials. These 
however are energetically higher and, unless we are in- 
terested in hybrid quarkonia with spin-exotic quantum 
numbers, not a dominant correction [83] to quark poten- 
tial model predictions. Obviously, hybrid meson mixing 
and decay is interesting in itself [84], and, in this con- 
text, a detailed study of the breaking of hybrid strings 
would be interesting. Finally, at light sea quark masses 
hadronic transitions between quarkonia, mediated by it 
radiation, become possible, the inclusion of which neces- 
sitates further modifications. 
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VII. SUMMARY 

We were able to resolve the string breaking problem in 
rif = 2 QCD, at one value of the lattice spacing aT 1 w 
2.37 GeV and of the sea quark mass, m < m s . 

To achieve this result, the systematic improvement of 
our methods beyond the latest lattice technology was cru- 
cial. In particular, we used highly optimized smearing 
functions to enhance the overlap of our test wave func- 
tions with the physical states, we employed an improved 
static action and we realized many off-axis source sepa- 
rations. We used stochastic estimator techniques (SET) 
to calculate all-to-all propagators. The variance of SET 
was reduced by exactly calculating the contribution from 
the lowest lying eigenmodes of 75 M, where M denotes 
the Wilson Dirac matrix. Further variance reduction 
was achieved by the (new) hopping parameter acceler- 
ation (HPA) technique. These methods, most of which 
are neither specific to the string breaking problem nor to 
the static charge sector of lattice QCD, can readily be 
applied to a large spectrum of problems. In general we 
would expect the gain factor from HPA to decrease at 
lighter sea quark masses while the improved convergence 
of TEA should compensate for this. 

We determined a mixing angle 0(r). Truncating the 
Fock basis after states containing four quark operators, 
the QQ component of the physical ground state was given 
by sin# and the BB content by cos9. We distinguished 
between the explicit and the implicit detection of mix- 
ing effects: a non-vanishing transition element between 
QQ and BB states is an explicit signal of mixing. Ad- 
ditionally, at r < r c the ground state energy will dom- 
inate, even within the BB operator, at large Euclidean 
times (implicit mixing). We were able to verify this be- 
haviour. At large r > r c , the lowest lying state will have 
a mass slightly smaller than twice the static-light mass, 
2m,B, and dominantly couple to the BB operator. We 
were unable to detect this signal in the QQ sector alone. 
Based on our mixing analysis, we expect such an im- 
plicit detection of string breaking to be almost hopeless, 
as high precision Wilson loop data at Euclidean times 
t = 0(5 fm) would be necessary — and this at distances 
r > 1.2 fm! 

We defined two string breaking distances: r c — 
15.00(8)a w 1.248(13) fm denotes the distance at which 
the energy gap between the two levels AE(r) = ^(r) — 
E\(r) assumes its minimal value AE C = 0.022(l)a~ 1 w 
51(3) McV while r s = 14.95(12)a w 1.244(16) fm denotes 
the distance of perfect mixing between the two states, in 
terms of the mixing angle 0(r s ) — 7r/4. The conversion 
into physical units has been made by setting ro = 0.5 fm. 
Note that r c — r s — 0.053(53)a is compatible with zero, 
within a standard deviation of less than 5 x 10 -3 fm. 

We would expect AE C w 51(3) MeV to increase with 



lighter sea quark masses and hence this value should be 
regarded as a lower limit to the case with massless or very 
light sea quarks. In real QCD also the string breaking dis- 
tance should decrease. We estimate r c = 1.13(10)(10) fm 
for this case. The first error is statistical and from the 
chiral extrapolation of previous results, the second in- 
corporates possible finite lattice spacing effects, the scale 
uncertainty in r = 0.5 fm and effects due to the inclu- 
sion of a third active sea quark flavour. The qualitative 
situation is depicted in Figure 22. 

We can define a transition rate g(r) = 
AE(r) sin[20(r)]/2 between QQ and BB states and 
note that g(r s ) = AE(r s )/2 and hence g(r c ) w AE c /2. 
In the lar ge N c l imit this means that at leading order, 
AE C oc yjrif /N c , if we are interested in the screening 
of a fundamental string by a sea of n/ massless flavours 
of fundamental particles, e.g. scalars or quarks. For the 
breaking of an adjoint string into two gluelumps the 
expectation reads, AE C cx 1/N C . In view of the precision 
of the n/ = 2 QCD results presented here, it should 
be worthwhile to dedicate renewed effort onto string 
breaking studies of SU(N C ) gauge theories with and 
without Higgs fields, to confirm the expectations, and to 
explore the applicability of large N c arguments to strong 
hadronic decays. 

We conclude that our study constitutes an important 
step towards the understanding of mixing effects and 
strong decays in quarkonium systems [78, 82]. Studying 
the energy between pairs of static-light mesons can also 
be viewed as a milestone with respect to a future calcula- 
tion of AqAq forces, which are related to nucleon-nucleon 
interactions [85]. 
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